function [me, rawspec] = gamma_measure(varargin)
%
% function to measure gamma oscillation (energy in 30-80 Hz band)
% from voltage traces.
%
% color scheme:
% raw signal: black or grey
% 30-80 Hz signal: blue
% 5-15 Hz signal: red
%
global VOLTAGE DFILE CONTROL
grey = [0.5, 0.5, 0.5];

testmode = 0; % if this is 1, we do a test
me=[];
rawspec = [];
% first, low-pass filter the VOLTAGE traces and then subsample them
nt0 = 0;
nt1 = 9999;

lpf = 250; % set lpf to 300 Hz

if(testmode == 1)
    newfigure('gammameasure', sprintf('Gamma Measure testdata'));
    % generate test data.
    % 10 seconds at 1 kHz sampling, with 2 secs noise, 2 secs at 10 Hz, and
    % 2 secs at 40 Hz, and 2 secs 10 + 40 Hz, 2 secs noise.
    rate = 50;
    nchan = 2;
    subsample = 10;
    fs = 1000000/(nchan*rate);
    fsr = fs/subsample;
    sV = [1, 100000];
    ratesub = rate * subsample;
    vo = -60 + 0.5*randn(sV(2), 1); % amplitude 0.5, at -60 mV
    sr = 1/fs; % sample rate, in time per point, per channel
    twos = floor(2*fs);
    ti = 0:twos-1;
    t = ti*sr;
    vo(ti+twos) = vo(ti+twos) + 1.5*cos(2*pi*10.*t)';
    vo(ti+2*twos) = vo(ti+2*twos) + 1.5*cos(2*pi*40.*t)';
    vo(ti+3*twos) = vo(ti+3*twos) + 1.5*cos(2*pi*10.*t)' + 1.5*cos(2*pi*40.*t)';
    vo = detrend(vo);
    nt1 = length(vo);
    nt0 = 0;
    % nt1 = (DFILE.nr_points-1)*DFILE.nr_channel(1)*DFILE.rate(1)/1000;
else
    rate = DFILE.rate;
    rate = rate(1);
    nchan = DFILE.nr_channel;
    nchan = nchan(1);
    subsample = 10;
        nt0 = 0;
        t=make_time(DFILE);
        nt1 = max(t(1,:));
    fs = 1000000/(nchan*rate); % convert to Hz.
    fsr = fs/subsample;
    %   sv = size(VOLTAGE);
    ratesub = rate*subsample;
    rec = varargin{1};
    sf=getmainselection;
    if(sf > 0)
        [s] = block_info(sf);
        if(strmatch(s.Name.v, 'mc-shock'))
            nt0 = 300;
            nt1 = 2000;
            CONTROL(sf).protocol = 'mc-shock'; % update it because it is probably wrong
        end;
        spike_thresh = CONTROL(sf).thresh;
    else
        spike_thresh = 0;
    end;
    it0 = floor(nt0/(rate*nchan/1000))+1;
    it1 = floor(nt1/(rate*nchan/1000))+1;
    if(any(VOLTAGE(rec, it0:it1) >spike_thresh))
        return;
    end;
    vo = pdetrend(t(rec, it0:it1), VOLTAGE(rec,it0:it1), 3);
    newfigure('gammameasure', sprintf('Gamma Measure on %s [%d-%d]', ...
        DFILE.filename, DFILE.frec, DFILE.lrec));
end;

[b,a] = butter(8,lpf/(fs/2));
Hd = dfilt.df2t(b,a);          % Direct-form II transposed structure
y = filter(Hd,vo);

% now subsample the data
t0 = t(1,:);
t1 = nt0:ratesub*nchan/1000:((nt1)*ratesub*nchan/1000);
y1 = interp1(t0, y, t1);


% bottom plot: signals and filtered versions
subplot('Position', [0.1, 0.1, 0.55, 0.25]);
plot(t0, vo, 'color', grey);
hold on;
plot(t1, y1, 'k.', 'markersize', 0.5);

% now calculate power in 30-80 Hz band
% make an arbitrary filter. Note f is specified in nyquist units...
nyf = fsr*0.5;
% f = [0 10/nyf 20/nyf 30/nyf 30/nyf 80/nyf 80/nyf  100/nyf 200/nyf 1];
% m = [0 0 0 0 1 1 0 0 0 0];
% [b, a] = yulewalk(32, f, m);
% Hd2 = dfilt.df2t(b,a);
% try chebyshev filter with min passband ripple
nord = 4;
wst = [30 80]/fsr*2;
[b, a] = butter(nord, wst);
%Hd2 = dfilt.df2t(b, a);
y2 = filter(b, a, pdetrend(t1, y1, 3));
plot(t1, y2, 'b-');

% low freq (5-15 Hz)

% f2 = [0 3/nyf 5/nyf 15/nyf 20/nyf 1];
% m2 = [0 0 1 1 0 0];
% [b2, a2] = yulewalk(32, f2, m2);
wst2 = [5 15]/fsr*2;
[b2, a2] = butter(nord,  wst2);
y3 = filter(b2, a2, pdetrend(t1, y1, 3));
plot(t1, y3, 'r-');

[h,w] = freqz(b,a,500);
[h2,w2] = freqz(b2,a2,500);
subplot('Position', [0.7 0.1, 0.25, 0.25]);
%semilogx(f*nyf,m,nyf*w/pi,abs(h),'r--');
semilogx(nyf*w/pi, abs(h), 'b--');
hold on;
%semilogx(f2*nyf,m2,nyf*w2/pi,abs(h2),'c--');
semilogx(nyf*w2/pi, abs(h2), 'r--');
set(gca, 'Xlim', [1 1000]);
%set(gca, 'Ylim', [0.0001*u, u]);
set(gca, 'XMinorTick', 'on');
set(gca, 'YMinorTick', 'on');
set(gca, 'Xtick', [1, 2, 3, 5, 10, 20, 30, 50, 100, 200, 300, 500, 1000]);
set(gca, 'Xticklabel', {'1', '', '3', '', '10', '', '30', '', '100', '', '300', '', '1000'});
grid on;
%
% envelope
yh1 = hilbert(y1);
Envelope1 = sqrt(imag(yh1).^2 + y1.^2);
yh2 = hilbert(y2); % compute hilbert on filtered data
Envelope2 = sqrt(imag(yh2).^2 + y2.^2);
yh3 = hilbert(y3); % compute hilbert on filtered data
Envelope3 = sqrt(imag(yh3).^2 + y3.^2);
Envelope0 = sqrt(Envelope1.^2 - sqrt(Envelope2.^2 + Envelope3.^2));
envref = (Envelope0 + Envelope2 + Envelope3 );
subplot('Position', [0.1, 0.4, 0.55, 0.25]);
plot(t1, Envelope1, 'color', grey);
hold on;
plot(t1, Envelope2, 'b-');
plot(t1, Envelope3, 'r-');
me.base = mean(Envelope0);
nEnvelope0 = Envelope0 ./ envref;
nEnvelope2 = Envelope2 ./ envref;
me.gamma = mean(nEnvelope2);
nEnvelope3 = Envelope3 ./ envref;
me.theta = mean(nEnvelope3);
subplot('Position', [0.1 0.7, 0.55, 0.25]);
plot(t1, nEnvelope0, 'color', grey);
hold on;
plot(t1, nEnvelope2, 'b-');
plot(t1, nEnvelope3, 'r-');
fprintf(1, 'Mean Envelope: %8.3f   5-15 Hz: %8.3f  30-80 Hz: %8.3f\n', ...
    me.base, me.theta, me.gamma);
f40 = length(find(nEnvelope2 > 0.5))/length(nEnvelope2);
f10 = length(find(nEnvelope3 > 0.5))/length(nEnvelope3);
fprintf(1, 'Fraction of time 10 Hz: %6.3f,  40 Hz: %6.3f\n', f10, f40);

%if(nargin > 0) % skip spectrum stuff if we are called in a loop
%    return;
%end;

% and finally, we need some raw spectral analysis of this data
Hs=spectrum.mtm;
subplot('Position', [0.7, 0.4, 0.25, 0.25]);
Hpsd1 = psd(Hs,pdetrend(t1, y1, 3),'Fs',fsr/1000);
Hpsd2 = psd(Hs,pdetrend(t1, y2, 3),'Fs',fsr/1000);
Hpsd3 = psd(Hs,pdetrend(t1, y3, 3),'Fs',fsr/1000);
max(Hpsd1.Data)
max(Hpsd2.Data)
max(Hpsd3.Data)
loglog(Hpsd1.Frequencies*1000, Hpsd1.Data, 'color', grey);
hold on
rawspec = Hpsd1;
loglog(Hpsd2.Frequencies*1000, Hpsd2.Data, 'b-');
hold on
loglog(Hpsd3.Frequencies*1000, Hpsd3.Data, 'r-');
if(nargin > 1)
    ref = varargin{2};
    if(~isempty(ref))
        loglog(ref.f*1000, ref.sp, 'k-');% averave reference (built up over time)
    end;
end;
set(gca, 'Xlim', [1 1000]);
%u = max([Hpsd1.Data; Hpsd2.Data; Hpsd3.Data]);
set(gca, 'Ylim', [1e-4 4]);
set(gca, 'XMinorTick', 'on');
% set(gca, 'YMinorTick', 'on');
set(gca, 'Xtick', [1, 2, 3, 5, 10, 20, 30, 50, 100, 200, 300, 500, 1000]);
set(gca, 'Xticklabel', {'1', '', '3', '', '10', '', '30', '', '100', '', '300', '', '1000'});
grid on;

function yd = pdetrend(x, y, n)
% polynomial subtraction for detrending. Order n on data x, y.
p = polyfit(x, y, n);
yd = polyval(p, x);
yd = y - yd;

